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We study the dynamics of parallel tempering simulations, also known as the replica exchange 
technique, which has become the method of choice for simulation of proteins and other complex 
systems. Recent results for the optimal choice of the control parameter discretization allow a 
treatment independent of the system in question. Analyzing mean first passage times across control 
parameter space, we find an expression for the optimal number of replicas in simulations covering 
a given temperature range. Our results suggest a particular protocol to optimize the number of 
replicas in actual simulations. 
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The effective simulation of proteins, glasses and sim- 
ilar complex systems has remained one of the defining 
challenges in computational physics. The main prob- 
lem in such simulations is slow relaxation due to bar- 
riers and bottlenecks. Parallel tempering - also known 
as the replica exchange method - promised a way out of 
this dilemma [H, 0, H- Here, canonical or generalized- 
ensemble simulations [J] are performed in parallel at dif- 
ferent values of a control parameter, most often the tem- 
perature. At certain times the current conformations of 
replicas at neighboring control parameter values are ex- 
changed according to a generalized Metropolis rule p|. 
An individual replica performs a random walk in control 
parameter space, allowing it to enter and escape local 
free energy minima. As a consequence, the state space 
is explored more evenly, especially e.g. at low tempera- 
tures. 

Replica exchange simulations are usually performed on 
massively parallel machines, with one or more comput- 
ing nodes dedicated to performing the simulation of a 
replica at a particular control parameter value. In order 
to optimize the use of resources, an important question 
is whether an optimal choice for the number of replicas, 
— or, equivalently, control parameter values [|| — exists, 
and how it might be determined. To our knowledge, no 
systematic investigation of this particular problem has 
been done to date. 

A likely reason is that this question is closely connected 
to the meta-dynamics of parallel tempering, and a full 
understanding of this method for complex systems - 
exhibiting broken ergodicity Q — is still missing [Tel . It 
is also not a well-posed problem. Apart from the num- 
ber of replicas, the main adjustable parameter is the 
distribution of control parameter values. The optimal 
number of replicas will depend strongly on the strategy 
used for the temperature discretization. Usually, a con- 
stant discretization is employed. However, often bottle- 
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necks exist in some control parameter regions, and speed- 
ing up the equilibration of the system is possible only by 
using a finer discretization in these regions. In order to 
investigate the question of the optimal number of repli- 
cas systematically, a method of discretization has to be 
employed that also complies with some optimality crite- 
ria. 

Major advances have been made recently in that direc- 
tion. Instead of concentrating on stationary distributions 
that arise from sampling, Trebst et al. [1, Q have investi- 
gated the flow across control parameter space and have 
provided an iterative scheme for adjusting the discretiza- 
tion to optimize the flow distribution. Subsequently, we 
have shown that optimizing the flow is equivalent to min- 
imizing the total first passage time to cross control pa- 
rameter space 10]. In this letter, we use those previous 
results as a basis for investigating the optimal number of 
replicas. We restrict ourselves to the situation of opti- 
mized flow and determine which number of control pa- 
rameter values — identical to the number of replicas — 
minimizes the first passage time of a single replica. 

In the following, we will consider parallel tempering 
with N+l replicas. Hence, we will assume N+l different 
control parameter values /3q < /?i < . . . < /3jv, i.e. we 
have N control parameter intervals [/3 n ,/3 n +i]. We will 
also use the conventions /3q — f3 m in and /3n — Pmax- For 
simplicity we will call the control parameter a (inverse) 
temperature in the rest of this paper. 

The time evolution of the probability P(n, t) that an 
individual replica is at temperature [3 n at time t can 
be approximated by a Master equation [ll[ in discrete 
time [l0| 

P(n, t + 1) = P(n,t) x 

[1 - W{fin -> /3 n -i) - W(/3 n /3„+i)] + 
P(n-l,i)W(fl,_i-.A l ) + 
P{n+l,t)W([3 n+1 ^f3 n ) , (1) 

where W((3 — > f3') are transition probabilities between 
neighboring temperatures. Of course, these probabilities 
depend on those temperatures, and the master equation 
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for replica exchange is characterized by symmetric tran- 
sition probabilities, 

W((3 (3') = W((3' ^(3) = W((3, /?') = W(f3', (3) . 

(2) 

The relation between stationary flow J and first passage 
time t in one-dimensional stochastic systems has been 
investigated in the context of channel flow in Ref. Il2l . 
Both quantities are related via 



J = C/t , 



(3) 



where C is a measure for the capacity of the channel. 
The mean first passage times for a single replica to cross 
the system defined by Eq. |lj in both directions is given 
by 0,0 
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while the channel capacity is simply the number of tem- 
perature values, C = N + 1 . It was shown in Ref. fiol that 

- for a particular number of control parameter values 

— the current is maximized (and therefore the first pas- 
sage time minimized) if the flow distribution is linear in 
the temperature number. This criterion allows an opti- 
mization of the temperature distribution ftj [lfj| . We will 
assume in the following that such an optimized distribu- 
tion of temperatures has been obtained. In this case, the 
effective transition probabilities in Eq. |T]) are constant 
across the chain of temperatures [Io| . 



W pt(/3i, A+i) = const for i = 0, 
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and this property will be essential in the analysis below. 

With increasing number of temperatures, i.e. finer 
discretization, the transition probabilities approach their 
maximum W — > Wq. Its numerical value depends on the 
particular implementation of the replica exchange algo- 
rithm and the choice of the exchange time scale. In this 
limit the mean first passage time shows the asymptotic 
behavior 
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i.e. it grows quadratically with number of replicas. How- 
ever, for smaller values of N the transition probabilities 
begin to decrease, leading eventually again to an increase 
in r for small N. We are interested in the value of N 
where r is minimal. 

This value will depend on the change of the transition 
probabilities with the control parameter interval, [/3,/3']. 
For the case of temperature as control parameter, this 
question has been investigated in depth 0, Hll [U [H, 
flQI ] . It has been found that the transition probability can 
be effectively approximated by [M E, EH, l2(il | 
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with f{x) a monotonically decreasing function obeying 
/(0) = 1. The important quantity here is b > 0, de- 
noting the scale of inverse temperatures over which the 
transition probability decreases. This scale is usually in- 
verse to the widths of the thermal equilibrium energy 
distributions at (3 and (3'. It decreases monotonically 
with system size and with the extensive heat capcity. In 
particular, it will be small near pha se transitions. More 
details can be found in Refs. [Mill, G3, GJ, OS 

In the following, we assume that for a particular system 
the functional form in Eq. (0 is the same over the full 
temperature range. The only dependence of W(/3j, ft+i) 
on the inverse temperature interval [/3j, /3j+i] is through 
the corresponding scale parameter that we denote by 

Under this assumption the requirement that - 
for the optimal temperature distribution — all effective 
transition probabilities are constant, Eq. ([5]), is equiva- 
lent to the condition that all individual arguments are 
identical. Hence, 
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f3% I = const = r 



holds. Introducing the average scale b 



N-l 



(8) 
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it can be seen that the constant r has the property 

rNb = (3 max - P min . (10) 

This property allows us to introduce the renormalizcd 
number of replicas: 



f3max (3n 



-N 



(11) 



which, in turn, allows us to cast the mean first passage 
time for a replica to cross the system into the parameter- 
free form 



(12) 



Minimizing v for a particular functional form f(x) of 
the transition probability decrease will give us finally the 
optimal number of replicas. 

We analyze Eq. (TT!?|) using the following functional 
forms for f(x), 



exp(— x) 
crfc {^Jn/Ax^j 
1-x 



(a) 
(b) 
(c) 



(13) 



In order to ensure they are comparable, we chose /(0) = 
1 = /'(0), i-e. the initial slope is identical for all three 
forms. Form (b) has actually been derived for tempera- 
ture intervals [13, [H, IS H3| and exhibits an exp(— x 2 )/x 
tail. It is the one most probable to occur in an actual 
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FIG. 1: Different functional forms for the decay of the tran- 
sition probability with control parameter distance, compare 
Eqs. © and (13]). 

situation. We have included forms (a) and (c) as worst 
case scenarios since they cover a wide range of behav- 
ior around (b): (a) exhibits a simple exponential de- 
crease, much slower than (b), while (c) exhibits a much 
faster, linear decrease; note that the latter is valid only 
for x < 1. Figure 1 shows a graphical comparison of the 
three functions. 

Minimizing Eq. (|12[) gives the optimal value for v, 

{ 1/2 : (a) 
v opt = { 1.05267 : (6) (14) 
I 3/2 : (c) 

These values are all of order one, despite the wide range of 
functional behavior of the transition probability decrease 
they describe. Rewriting Eq. (jlip we obtain our final 
result for the optimal number of replicas 

N op t — Vopt (Pmax ~ flmin) jb . (15) 

The ratio of the full temperature range of the simulation 
to the average scale b is the main determining quantity in 
that equation. In particular, it controls the order of mag- 
nitude for N opt . Particular functional forms for the dis- 
tance dependence of the transition probabilities appear 
to have less influence since v opt is of 0(1) only. Although 
the range of values given in Eq. (|14p for v opt is still cov- 
ered by factor of three, this is a worst case scenario, and 
actual values for realistic functional forms will be closer 

to (33b). 

The influence of the different functions, Eq. (fl3|) . is 
more important for the form of the minimum of the mean 
first passage time. Figure 2 shows how the mean first pas- 
sage time changes when the number of replicas N devi- 
ates from the optimal value N opt . In order to enable com- 
parison we use the renormalized replica number, Eq. (jlip . 
as variable here. While the minimum is pronounced for 
all functional forms, it is steepest for form (c), i.e. the 



fastest decreasing probability function, and most shallow 
for form (a). 

Before we discuss the consequences of Eq. (fT5"|) for pro- 
tocols to optimize the number of replicas, we need to 
address some subtleties of the above derivation that we 
skipped over in favor of a compact derivation: 

(i) Equation (TJ) is an effective description of the long- 
time properties of the random walk of replicas in par- 
allel tempering simulations. The transition probabilities 
for such a long-time description differ from observed ac- 
ceptance rates, which are, on the replica exchange time 
scale, short time properties. In particular, it has been ob- 
served that for optimized flow the observed acceptance 
rates are not constant 0] . This is due to broken ergodic- 
ity 0] at particular control parameter values, that gives 
rise to a hierarchical, tree-like structure for the random 
walk of replicas 10]. Observed flow and acceptance rates 
are just projections onto the Id chain of control param- 
eter values of the more complicated flow processes on 
the tree. However, the possibility of flow optimization 
shows that for the long-time transition probabilities the 
property (O holds nevertheless. Such discrepancies be- 
tween short-time and lon g-time p roperties of stochastic 
processes are well-known [lj, [l5|]. By using Eq. (7]) we 
implicitly assume that those effective transition probabil- 
ities exhibit the same qualitative behavior with control 
parameter difference as it was derived for the short-time 
transition probabilities. We feel that this is justified since 
our qualitative results for N opt , Eq. (|15p . are independent 
of the particular functional form. 

(ii) J and r differ by a factor of N + 1, see Eq. (3]). The 
mean first passage time is usually a good estimate for 
the lowest eigenvalue of the equation system ffl, i.e. it 
determines the time scale of equilibration [14], EH > Since 
we are interested in fast equilibration, r is the more ade- 
quate quantity than J to use for comparing systems with 
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FIG. 2: Dependence of relative mean first passage time, 
Eq (|12p on deviations from the optimal number of replicas 
for the different functions of Eq. (|13p . 
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different numbers of replicas and to optimize with respect 
to N. 

(iii) We have omitted a discussion of the replica exchange 
time scale. Depending on the frequency of replica ex- 
change moves the time scale of Eq. ([1]) may differ by 
a factor of N. However, our result (|15|) is stable with 
respect to changes in the power of N. Although the 
particular numerical values change, the value of v opt 
remains a 0(1) constant if the exponent in Eq. (TT2|) 
changes from two to one (v opt = 1,1.6671,2) or three 
(y ovt = 1/3,0.820024,4/3). 

What are the consequences of our results, particularly 
of Eq. (fT5|) , for protocols to optimize the number of repli- 
cas? The main result of Eq. (fl~5|) is that it identifies, sep- 
arates, and quantifies the contribution of various proper- 
ties of the simulation system to the optimal number of 
replicas. In particular, it exhibits the quantitative hier- 
archy of the individual contributions. The above analysis 
also shows the importance of how transition probabilities 
change with the control parameter interval. To our sur- 
prise, for the inverse temperature as control parameter, 
this complex contribution could be summarized formally 
into the averaged scale 6, Eq. Taking into account 
the dependence of the scales on the extensive prop- 

erties of a system [l(], [HI suggests that N opt scales 
with system size V as N opt cx \fV . 

We note that the determination of b in actual simula- 
tions is by no means simple. Since it is defined for the sit- 
uation of optimal control parameter spacing for a partic- 
ular number of replicas, such an optimization would have 



to be performed beforehand. Also, since it describes the 
behavior of the effective transition probabilities, see the 
above discussion, it would have to be determined from 
the flow distribution together with the actual first pas- 
sage time, upon slightly varying the discretization. 

Instead, our analysis suggests that the direct approach 
to optimize the number of replicas is the most promising 
one. Comparing first passage times of replicas to cross 
the simulation system for the optimized discretization is 
readily possible for different values of N. Figure 2, in 
particular form (b), can then be used as a guideline to 
extrapolate to N opt . 

In summary, we have studied the dynamics of paral- 
lel tempering simulations. Analyzing these dynamics, we 
have determined the main factors influencing the optimal 
number of replicas in such simulations and their quantita- 
tive hierarchy. Since the evaluation of the essential term 
6, the average scale of transition probability decrease, 
may need costly computations, we propose to base opti- 
mizing the number of replicas on the generic behavior of 
a replica's first passage time to cross the simulation sys- 
tem given in Fig. 2. The technique of replica exchange 
has become the method of choice for the simulation of 
proteins and other complex systems. The above results 
add to its understanding, and we believe they will also 
advance its practical use. 
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